4T1 Mets RNA Seq - Enrichment Analysis

library(Biobase)
library(SummarizedExperiment)
library(openxlsx)
library(stringr)
library(tidyverse)
library(ggplot2)
library(reactable)
library(hypeR)

PATH <- file.path(Sys.getenv("MLAB"),"projects/brcameta/exosome/4t1_brca/")
DPATH <- file.path(Sys.getenv("CBM"),"otherStudies/RNAseq/2023-03-22-YuhanExosomeBrCa")

do_save <- TRUE

Loading Signatures and Data

all_sigs <- readRDS(file.path(PATH,"data/signatures_symbol.rds"))
all_sigs <- all_sigs[lapply(all_sigs, length) != 0]

deseq_list <- readRDS(file.path(PATH,"data/deseq_list.rds"))

# PAMM Genelists
pamm_paths <- Sys.glob(file.path(PATH, "data/PAMM_Genelists/*"))
pamm_names <- lapply(pamm_paths, stringr::str_match, pattern="genelist-\\s*(.*?)\\s*.xlsx")
pamm_names <- unlist(lapply(pamm_names, function(x) x[2]))

pamm_genelists <- lapply(pamm_paths, openxlsx::read.xlsx, colNames = FALSE)
pamm_genelists <- lapply(pamm_genelists, function(x) x$X2)
names(pamm_genelists) <- pamm_names

hallmark_genesets <- msigdb_download(species="Mus musculus", category="H")

dat <- readRDS(file.path(PATH,"data/4T1_mets_sExp.rds"))
background <- unique(unlist(SummarizedExperiment::rowData(dat)$mgi_symbol))

#Other genesets
methyltransferase_genes <- read.delim(file.path(PATH, "data/methyltransferase_genes_ncbi.txt"))
dnmt_genes <- methyltransferase_genes$Symbol[methyltransferase_genes$Symbol %in% background]

hallmark_genesets$DNMT <- dnmt_genes
h3_meth_genes <- c("PRC1", "PRC2", "EZH2", "EED", "SUZ12", "RING1", "RNF2", "BRD2", "BRD3", "BRD4", "KDM1A", "KDM2A", "KDM3A", "KDM4A", "KDM5A")
library(babelgene)
mouse_h3_meth_genes <- babelgene::orthologs(genes = h3_meth_genes, species = "mouse")
hallmark_genesets$H3METH <- mouse_h3_meth_genes$symbol

hypeR enrichment

Hypergeometric Hallmark

hall_hyp <- hypeR::hypeR(all_sigs,genesets=hallmark_genesets,background=background, test="hypergeometric", fdr=0.05)
hypeR::hyp_dots(hall_hyp,merge=TRUE,top=15)

hypeR::rctbl_build(hall_hyp)

PAMM

pamm_hyp <- hypeR::hypeR(all_sigs,genesets=pamm_genelists,background=background,fdr=0.05)
hypeR::hyp_dots(pamm_hyp,merge=TRUE,top=25)

hypeR::rctbl_build(pamm_hyp)

if (do_save) {
  hypeR::hyp_to_excel(hall_hyp,file.path(PATH,"data/hallmark_hyp.xlsx"))
}

KS-based enrichment

## Add gene symbols to DESeq tables

deseq_list <- lapply(deseq_list, function(Z) {
  as.data.frame(Z) %>% 
    tibble::rownames_to_column(var="ensembl_id") %>%
    dplyr::mutate(geneID=unlist(rowData(dat)[ensembl_id,"mgi_symbol"]))
  })

## rank by up-regulation
rank_signatures_up <- lapply(
  deseq_list, function(sig) sig %>% 
    dplyr::filter(geneID!="") %>%
    dplyr::arrange(desc(log2FoldChange)) %>%
    dplyr::select(geneID,log2FoldChange) %>%
    tibble::deframe())
names(rank_signatures_up) <- paste(names(rank_signatures_up),"up",sep="_")

## rank by down-regulation
rank_signatures_dn <- lapply(
  deseq_list,function(sig) sig %>% 
    dplyr::filter(geneID!="") %>%
    dplyr::arrange(log2FoldChange) %>%
    dplyr::select(geneID,log2FoldChange) %>%
    tibble::deframe())
names(rank_signatures_dn) <- paste(names(rank_signatures_dn),"dn",sep="_")

rank_signatures <- c(rank_signatures_up,rank_signatures_dn)
rank_signatures <- rank_signatures[order(names(rank_signatures),decreasing = TRUE)]

Hallmarks + Pamm

all_genesets <- c(hallmark_genesets, pamm_genelists)
max_fdr <- 0.05
ks_hall <- hypeR::hypeR(signature=rank_signatures,genesets=all_genesets,test="ks",fdr=max_fdr,plotting=TRUE)
hyp_dots(ks_hall,merge=TRUE,fdr=max_fdr/5,top=25) + ggtitle(paste("FDR ≤", max_fdr/5))

hypeR::rctbl_build(ks_hall, show_hmaps=FALSE)

IS C UP

print(ks_hall$data$is_c_up$plots[ks_hall$data$is_c_up$data$fdr<=0.01])

## $HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION

## 
## $HALLMARK_INTERFERON_GAMMA_RESPONSE

## 
## $HALLMARK_KRAS_SIGNALING_DN

## 
## $HALLMARK_XENOBIOTIC_METABOLISM

## 
## $HALLMARK_TNFA_SIGNALING_VIA_NFKB

## 
## $HALLMARK_BILE_ACID_METABOLISM

## 
## $HALLMARK_INFLAMMATORY_RESPONSE

## 
## $HALLMARK_OXIDATIVE_PHOSPHORYLATION

## 
## $HALLMARK_MYOGENESIS

IS C Down

print(ks_hall$data$is_c_dn$plots[ks_hall$data$is_c_dn$data$fdr<=0.01])

## $HALLMARK_MYC_TARGETS_V1

## 
## $HALLMARK_E2F_TARGETS

## 
## $HALLMARK_G2M_CHECKPOINT

## 
## $HALLMARK_CHOLESTEROL_HOMEOSTASIS

## 
## $HALLMARK_MYC_TARGETS_V2

## 
## $HALLMARK_MTORC1_SIGNALING

## 
## $HALLMARK_MITOTIC_SPINDLE

## 
## $HALLMARK_UNFOLDED_PROTEIN_RESPONSE

## 
## $`Mouse Cell Cycle`

IR IS Up

print(ks_hall$data$ir_is_up$plots[ks_hall$data$ir_is_up$data$fdr<=0.01])

## $HALLMARK_G2M_CHECKPOINT

## 
## $HALLMARK_E2F_TARGETS

## 
## $HALLMARK_MYC_TARGETS_V1

## 
## $HALLMARK_MITOTIC_SPINDLE

IR IS Down

print(ks_hall$data$ir_is_dn$plots[ks_hall$data$ir_is_dn$data$fdr<=0.01])

## $HALLMARK_MYOGENESIS

## 
## $HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION

## 
## $HALLMARK_OXIDATIVE_PHOSPHORYLATION

## 
## $HALLMARK_INTERFERON_GAMMA_RESPONSE

## 
## $HALLMARK_ADIPOGENESIS

IR C Up

print(ks_hall$data$ir_c_up$plots[ks_hall$data$ir_c_up$data$fdr<=0.01])

## $HALLMARK_XENOBIOTIC_METABOLISM

## 
## $HALLMARK_TNFA_SIGNALING_VIA_NFKB

## 
## $HALLMARK_INTERFERON_GAMMA_RESPONSE

IR C Down

print(ks_hall$data$ir_c_dn$plots[ks_hall$data$ir_c_dn$data$fdr<=0.01])

## $HALLMARK_MYC_TARGETS_V1

## 
## $HALLMARK_MYC_TARGETS_V2

## 
## $HALLMARK_MTORC1_SIGNALING

## 
## $HALLMARK_CHOLESTEROL_HOMEOSTASIS

## 
## $HALLMARK_E2F_TARGETS

## 
## $HALLMARK_UNFOLDED_PROTEIN_RESPONSE

## 
## $HALLMARK_G2M_CHECKPOINT

## 
## $HALLMARK_MYOGENESIS

## 
## $HALLMARK_TNFA_SIGNALING_VIA_NFKB